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ABSTRACT 

A method is developed, using the quadratic form of canonical un- 
coupled state variables as a Lyapunov function to determine the switch- 
ing logic for a quasi-optimum control of a non-linear dynamic systema 
Analytic arguments are presented to support the method, which is capable 
of being extended to any order system and any number of controls, sub- 
jeot to certain practical limitations noted in the analysis. 

A computational scheme for determining the canonical state variables 
and control functions is presented and a topological interpretation is 
made of the choice of variables used for the control logic calculation. 

The controller based on the method described is applied to a none 
linear second-order system. Phase trajectories for both the uncontrolled 
and controlled systems are obtained by means of a digital computer simula- 
tion. Various aspects of the theoretical limitations of the method pre-e 
sented are investigated and the experimental results are analyzed with 
respect to the theoretical predictions. 


#The system is postulated to be described by a system of differential 
equations of known forme 
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1. Introduction. 

The stability of dynamic systems whose equations of motion are such 
that their solutions are difficult or impossible to obtain in closed form 
may successfully be analyzed by means of Lyapunov’s Second Methed. 

We have a given dynamic system described by the set of differential 
equations as 

x ja g(x); g(0) = 0. (1-1) 
Then, referring to Figo (l-l), the origin, according to LaSalle and 
Lefschetz” is: 

stable whenever for each RCA there is an r ¢R such 

that if a path (a motion) g~* initiates at a point x° of 

the spherical region S(r) then it remains in the spherical 

region S(R) ever after; that is, a path starting in S(r) 

never reaches the boundry sphere H(R) of S(R)3 

asymptotically stable whenever it is stable and in addition 

every path g* starting inside some S(R,), Ry» o, tends to the 

origin as time increases indefinitely; 

unstable whenever for some R and any r, no matter how small, 

there is always in the spherical region S(r) a point x 


ae) that the path e* through x reaches the boundary sphere 
H(R). 





Fig. (1-1)° Geometric Interpretation of the Various Types 
of Stability. 








Now suppose we have a given dynamic system described by the set of 


differential equations written in vector form as 


x = g(x,u) (1-2) 

xX = Fx + Du + g (x,u) (=3) 
where F is a constant square matrix, D is the control distribution 
matrix, u is the control function vector, and e (x,u) contains terms of 
a second and higher power in x and ue The linearized form of the equae 
tion is 

x = Fe + Du (1-4) 
and can be said to adequately describe the motion of the system for 


small values of Xe 


According to Kalman and Bertram“, a Lyapunov function may be de- 


fined as some scalar function of the state variables V(x) with the fol- 


lowing properties: (a) V(x) )0, V(x) ¢ O when x + X,, where is the 


<% 
equilibrium state, and (b) V(x) = V(x) = Owhenx = X,0 

We wish to apply Lyapunov's Second Method to the design of a cone 
troller so as to return the state vector (x) to a position of equilib- 
rium (x,) in some ontimum manner. No loss of generality will result 
if we choose (x2) to be the origin of the state space, and henceforth 
the equilibrium position will be considered as such. 

Lyapunov stated that if a Lyapunov function V(x) exists (being 
positive definite) in some neighborhood of the origin, and if ~V(x) 
is also positive definite in some neighborhood of the origin, the 
system is asymptotically stable. 

In general the method employed will be to make V(x) as negative as 


possible in order to drive the system to the origin as quickly as 


possible. 








The Second Method of Lyapunov leaves open to the individual inves- 
tigator the choice of a particular Lyapunov function to be used in the 
analysis of a given system. For a given system then, there are possibly 
an unlimited number of such functions. The quadratic form of the canoni- 
cal state variables as the particular Lyapunov function employed in 
the design of the controller was governed by the simplification obtained 
in the control logic. 

Analysis and design for a controller of a linear, stationary system 
using the quadratic form procedes in a straightforward manner. In the 
case of a non-linear system however, the limitation that the linear ap- 
proximation is opie eloficnteily accurate only within some region surrounding 
the equilibrium position poses a rather severe restriction on the extent 
of the state space for which a given control can be useful. 

In order to overcome this limitation to some degree the controller 
was designed in such a manner that the control logic was obtained by 
successive linear approximations to the actual system at points on the 
system trajectory. 

For best performance of the system the time involved in the calcula- 
tion of the control logic should be infinitesimal so that there would be 
no appreciable time delay in the application of the proper control func- 
tion. That is, the system state vector would not have moved from the 
point (x, ), where the state variables were sampled, to a new position on 
the trajectory, (xj) where the control function was actually applied. 

In actual practice a finite euniguiaeion time is of course necessary 
so that the state vector will always have moved from the point of sampling 
before the new value of control function can be applied. 

For plants with long time constants, such as might exist in chemical 
or metallurgical processes or ship stabilization systems, presently 


%) 





available digital computers have computation speeds adequate to success= 
fully provide on line, real time control. 

Application of this method to plants with short time constants will 
depend upon increased speeds of the computers and the development of more 


refined methods of calculation. 








Lie Theory. 


If we have given a Lyapunov function V(x) = x! Px we may assert 
that* the equilibrium state x, 2 O of a linear dynamic system 
Xs Fx . (2-1) 
is asymptotically stable if and only if given any symmetric positive-= 
definite matrix Q, there exists a symmetric, positive definite matrix P 
which is the unique solution of the set of linear equations described by 
Flip 4 PF = -Q. (222) 
A simple procedure for calculating the control logic attributed to 
R.W. Bass? procedes as release: Choose Q>0 arbitrarily. P>O may then 
be calculated from (2-2). ||x|/@P will then be a Lyapunov function for the 
system (2-1). We now choose u(t) so as to make V (x) as negative as 
possible. Let 
Vix) = xi px = \Ix||@P. then 
V(x) se xnPx 4) Xora (2-3) 
We wish to consider the effect of control on V. Taking the transpose 


of (1-4) we obtain 


g°oe xlpT 4 ulpT (224) 
and 
Vix) = xi pl px + x! PFx + ulplpx + xlPpu 
or 


V(x) 


xi (Flip + PF)x 4+ 2u'p Px ° 


PD. Simplifying this gives 


Since P is symmetric, pip 


V(x) = -xlqQx + 2ulpTpx 
or V(x) = -[x]/Q + 2utdTpx . (2-5) 


To make V(x) as negative as possible we want the term 2u'D! Px to be 


as negative as possible. Therefore the control logic is given as 








ui(t) = -ay sgn (DI Px); (2-6) 


where the a: 


; are constants denoting the maximum allowable magnitudes of 


each individual control effort uj. It is obvious from (2-6) that in 
order to make the term 2u'pl Px as negative as possible the individual 
controller u,(t) must always exert a maximum effort, and in such a 
direction that the term 
“a; sen (p'Px), : 

is, in fact, negative. Thus the control logic becomes merely logic for 
detecting the appropriate switching points of the "bang-bang" controller 
a:(t). 

In the case where there is only one controller the vector u may be 
of the form such that ae wy (OpOSO:¥% ws. u). The system equation is 


then of the form 


0 1 om, % « 90 0 
0 0 1 0 
x = - x + : (a7) 
1 ° 
ae -b, | u 


where the b; are the coefficients of the characteristic equation of the 
unforced system, taken in reverse order. 
If we make the equiaiienee transformation Pron the physical state 

Space to a canonical state space denoted by the matrix equation 

y= & (2-8) 
then by differentiating we get 

y 2 Gt (29) 
and equation (1-4) becomes 


1 


y =..¢ Fey + -1Du . 








We choose G such that 
Ging = A. (210) 
where AL is a diagonal matrix with the same eigenvalues as Fi 
For the system with one controller such as described by (2-7) we 
may further choose the particular transformation matrix G@ such that 
w’Du = B 
where EE = (a, a, 8, .. 08). This has the effect of weighting the 
value of all the state variables equally in determining the sign of u(t). 
If the P of (2-2) is chosen to be the identity matrix I and the AW 
matrix substituted for F, equation (2-2) becomes 
AlT1 + Le =o. (2-11) 
If all the real parts of the eigenvalues of F are negative, so then are 
all those of /V_, and Q is positive definite. Also 
Wiz) = gCcA tT on TA. )y 4 opty. (2-12) 
The first term on the right side of (2-12) is again negative and the 


control logic is given by the second term. For a single controller, 


‘we have 


u(t) = -a sgn ( ¥ Ys) (2215) 
In other words the switching function is such that the control must 
change sign whenever the sum of the n weenie? state variables changes 
sign. 


Consider the given form 
ys AY + @'Du. 


#It is always possible to perform the equivalence transformation (2-10) 
provided F has no repeated eigenvalues. A system described by an nth order 
differential equation having repeated roots may be transformed into a pare 
tially uncoupled set of lst order differential equations. It will be 
assumed that the systems considered in this paper have no multiple roots. 
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If we take the LaPlace transform we get 

sY(s) - y(0) = A.Y(s) + G-1Du(s) 
or 

YGs) (el = FS) we yO) an G--Duge (214) 
Since from (2-6) the control variable should assume only values of + a 
or -a , we stipulate that it is a constant during a given control 
period (between switchings). We then have the laPlace identity 

(G-1Du)(s) = =. 

S 


Solving (2-14) for a given component y, we get 


Yi(s)(s -2,) = yy(0) + E 





or 
yi (0) a5 


——— ae =-15 
(s -d,) x s(s -),) tog 


Y3(s) 


Then the time domain solution of (2-15) is 











yy(t) = ys(o)e*#® - 2 (1 -9*8") 
a 
or 
a of ri% a 
y,(t) = fy3(0) + er be : — (2-16) 


An interesting result becomes apparent from this solution as t ine 
creases without bound. For the case where A4 is negative (a consequence 


of Q being positive definite) we see that 


Cee: a. 
y,(t) => ap Cima t+ -» ©0 


This result indicates that although this form of controller is useful in 
making V(x) more negative in an asymptotically stable system, it will not 
result in V(x) >0O as t—~ oO if u is a function of the type seen in 


Mig. (2-1). 





= — = se 








Fig. (2-1) Control Function U, As An Ideal Relay. 


In fact it will result in a condition of chatter-motion around the origin 
as the controller becomes dominant at points in the state space where the 
Yi <|lel]: where |le|| assumes some particular small value. 

This objectionable feature may be overcome in the case of asymptoti- 
cally stable systems by choosing the control function of the form shown 
in Fig. (2-2) below so that the controller is effectively disengaged 
within some region arbitrarily close to the origin, allowing the system 


to reach the equilibrium condition in an unforced manner. 





Fig. (2-2) Control Function U; As A Relay With Dead Zone. 


If the uncontrolled system were such that one (or more) roots con- 
3 


tained positive real parts, (in this case =-Q can no longer be positive 


definite) then for some particular points in state space, it is still 








possible that 
° 2 
V(x) = -||x llQ ° 
However V(x)>0 and remains SOc 


Consider a plant governed by the equations 
dx/dt = Fx 4 Du(t) 
where the control variables are subject to the constraints, 
luz (t)] Side < OO (i = 1, som symm). Gee27) 
This is essentially the relay or saturating servo problem. 
The problem is to return every initial state to the origin. 
It is well known that if F has sigenvalues with posi- 
tive real parts then there are some states which cannot be 
returned to the origin by any control subject to the cone 
straints (2-17)*. 


If a controller of the form described is applied to the unstable 
system, it is apparent from (2-5) that control may be maintained so long 
as the inequality (2-18) below exists.* That is 

-|x]#2 + 2utptpx < o. (2-18) 
From the solutions (2-16) the behavior of each y;(t) for 0<t<k, 


where k denotes the switching period is as in Fig. (2-3). 


¥, (t) 





Fig. (2-3) Plot of y;(t) vs. Time For u, = <a; sgn (By), 
Stable Eigenvalue. 
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Figs (2-4) Plot of y;(t) vs. Time For u; s -a; sgn (Gana, 
Unstable Eigenvalue, y;(0)>| i 
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| y, (t) 


4 "i (0) 





Figs (2-5) Plt of yi (t) vs. Time For u; = -a; sgn Cao 
Unstable Eigenvalue, yil0)<| ay | . 
Ai 
In the event that the system distribution matrix D were such that 
the controls were uncoupled, i.e. for each yy there were a corresponding 
u;, then it would be possible to design a controller whose individual 
elements changed sign at the zero crossings of the individual state 


variables. 


uk 








From Fig. (2-5) it can be seen that the element would switch at 
t = T. It is possible to calculate T, = T,y(yq(0)) from (2-16) for 


the case where 4,50, y,(0)>0, u, <0 » 








re 
y,(t) * {¥1(0 “ an _— * = » let y,(t) = 0 














Ad i 
then 
a. ‘XT a, 
{3i(0) - i wo eee 
dey Des 
giving 
a 2 P a { dy 
Ba U Aayy(0)-a4 
and 
(2-19) 
T = ame 2 eee 2-19 
dy ’ { eA et 


From (2-19) it is clear that for y,;(0) =» 0, T =& —~ in(l) = O, 
indicating that this system also will result in shat Coecilaelian around 
the origin. 

Selection of a control function uz of the form of Fige (2-2) would 
be unsatisfactory in this case because the uncontrolled system is une 


stable at the origin. A better form of function would be that of the 


form of a saturating emplifier with a high gain. (See Fig. (26)). 


Fig. (2-6) Control Function u; With the Characteristics of a 
Saturating Amplifier. , 
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Employing LaSalle and Lefschetz's definitions® 


of stability des 
cribed in Section 1., a system whose eigenvalues all have negative real 
parts is asymptotically stable with no control, merely stable with a 
"bang-bang" controller (Fig. (2-1)), and again asymptotically stable if 
the controller contains a dead zone as in Fig. (2-2). Given an arbit- 
rary H(R), a system with eigenvalues lying in the right half plane is 
defined as unstable, and the same system with a controller of the form 
described may be defined as stable if X° is restricted to lie within 
some S(r) whose redius r is of some necessarily small value. The 
restrictions on the allowable values of r are necessarily related to 
‘the constraints on the magnitudes of the control functions (2-17). 
Given a non-linear system described by the equation x & F(x)X5 
the linear approximation of the non-linear system referred to by Kalman 


and Bertram® may be obtained by a Taylor series expansion of the coef= 


ficient matrix F(x) at the origin. The linearized, constant, matrix F 
OF; (x) ; 
' o Xi 
with respect to x at x = x(0) = OQ). 


has the elements (The matrix F is the Jacobian of F(x) 


Suppose a system described by the equation 

% = F(x)x + Du 

.is to be controlled in an envircmment where the disturbances are of 
such a magnitude as to carry the state variable beyond the region of 
validity of the linear approximation given by (1-4). Referring to 
Fig. (2-7), this condition corresponds to an x° lying outside of S(r). 
By expanding F(x) about the point x° and obtaining the Jacobian of 
F(x) with respect to x at x = x° we may obtain a linear approxima- 


tion at x° 


(2.20) 
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Figo (2-7) Geometric Interpretation of Successive Linear 
Approximation Process. 

Control logic based on this approximation will be valid while the 
trajectory remains within the spherical region S(p)o If S(p") is the 
region of validity of the linear approximation due to the nth sempling » 
then by an iterative process such that the trajectory from X27 to xn+t 
always remains inside S(p™) we may obtain an optimum control policy 
based on the particular Lyapunov function employed for the complete 


trajectory for x° tox, = 0. 
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3. Application of Theory to Design of a Single Controller. 

It was shown in Section 2. that the control logic resulting from 
the choice of the quadratic form of the uncoupled state variables as 
the system Lyapunov function is of the form 

u(t) = -a sen} 2, vit 
when the system has only one controller. Thus, the control problem is 
reduced to determining the sign of the sum of canonical state variables 
and insuring that the control, u(t), assumes the opposite sign. The 
steps involved in this process consist of the following: 

1. Measuring the physical state variables (x). 


2 Calculating the canonical state variables (y) from the 
measured (x). 


Se Forming > Yie 
4. Applying the proper control u(t) = -<=a Sgn) Vie 
As was shown in Section 2., the variables (y) may be calculated by 
the equation 
y # gntx ° (3=1) 
The problem then reduces to evaluating Gang 
It may be shown that a matrix G such that 
Gee Sali Gs.) 
where _/\is a diagonal matrix, may be formed by calculating the eigen- 
vectors of the matrix F and constructing an array in which each eigen- 
vector is a particular column of the array. Once a G has been obtained 
it is merely necessary to calculate the inverse, q-l, in order to be able 
to calculate the state variables (y). 
The elements of the diagonal matrix + are the eigenvalues of F. 
If the system is oscillatory the elements of ALwill necessarily be 
complex. Since the elements of the coefficient matrix F of the original 


15 








differential equation are all real, consideration of (3-2) reveals that 
for an oscillatory system there must be complex slements in G@ and ule 

When Gis complex, care must be taken in evaluating Gel. If we 
let G= RG 4 JWG, and attempt to form aes = (@e)-t y (a G)~+ 
difficulty arises immediately. Since complex roots arise in conjugate 
pairs, there will be two eigenvectors (columns) of ((RG) due to the 
identical real parts of the conjugate pairs. These vectors will be 
identical, and therefore linearly related so that (GG) will be singular, 
Similar reasoning shows that (3G) will also be singular, so that it is 
impossible to obtain (®G)-1 and (g G)-i, 

Let us form a new matrix H in the following manner. If the order 
of G is n, the order of H will be 2n so that for every complex element 
of G there will be four elements of H arranged in a square pattern. The 
two elements on the principal diagonal of the squares have a value equal 
to the real part of the corresponding complex element of G. The ab- 
solute values of the remaining two elements are equal to the imaginary 
part of the corresponding complex element cf G, the element in the lower 
left hand corner, taking the positive sign, while that in the upper right 
hand corner, the negative sign. Where the corresponding element in G 
is real, zeros appear in the locations in H corresponding to the missing 
imaginary parts. 

If now we form H-!, then we may obtain the elements of G-1 from 
the proper locations of Hel, 

However since we are interested in obtaining the (y) corresponding 
to the measured (x), the simplest procedure is to form the 2 x n sup- 
plementary vector tion has alternately the real and imaginary 


parts of the physical state variables. Since these state variables may 


Ls 








be measured and have physical significance, they may only be real, so 
that a consists of the values of the state variables alternating with 
zeroSse Matrix multiplication of q-l on the right by =" results ina 
2 xX n column vector whose elements are the real and imaginary parts of 
the canonical space variables (y). 

EXAMPLE. The product of a complex 2 x 2 square matrix and a 
column vector with real elements 

Yi) | Pin omtyy Beat Stig — 


x 
Roy + Joy Roe + Ilo Xo 


y 


do 


Then 
Ba ley | ®11%1 ote nace een ed laze 
Yo Roi% + J2g)%1 + Rop% + Jlgo%p 
Rearranging gives 
¥y = (Rye + Rypxe) + §(Tqy% + Typ%q) Ce) 
Yo = (Ropx%, + Ropxe) + J(Ip)%, 4 Iggxq) (3-4) 


Suppose instead we forma 4 x 4 matrix of the form H7!, Then 


RY} Ruy tT. Fag = 51a 2 xy 
AY} | t32n. Far +5li2 Re aie? 
Ryz | Roy 5g Reg «= -JI 22 X2 
dyo}- | dle. Re, J1g2 Rep 0 
and = 
Ry, > Ry% + Ror, 
| Sy. = 3t1%1 + ST ygxg 


QYo = Rox, + Rox 


X¥2g = J1g 1% 


+ 


JlooxXo ° 
Rearranging gives 


Yy = (RQyXy + Ryoxe) + jlTy%, + = Ty2%) 
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ya = (Ro1%, + Roexe) + jllg}%, +  Tp9%e) 

which is identical to (3-3) and (3-4). 

It is evident from the above discussion that in the case of an oscile- 
latory system the transformation 

y = oly (3=5) 

may map a given physical state variable x; from a point on the real 
axis to some point in the complex plane. Thus for a given coordinate 
axis in the physical n-space there is a complex plane in the canonical 
SpACCe 

Since, however, there is a one for one correspondence between points 
in the two spaces, a control which succeeds in reducing all state vari- 
ables to zero in the canonical space will have also reduced the real 
state variables to zeroe 

The computation method employed to determine the transformation 
matrix G was based on a computer program MATSUB, programmed by 
Louis W. Ehrlich of the University of Texas, which calculates eigen- 
values and eigenvectors of a matrix up to order 50. The general method 
employed is to make an initial guess for each eigenvector of the matrix 
F and to converge on the correct value by an iterative scheme due to 
E.E. Osborne®, The original program employed a guess of 1.0 for all com- 
ponents of the eigenvectors. Since it was assumed that the system 
matrix F was not changing at more than a moderate rate, a modification 
was made to the program so that the most recently calculated value for 
each eigenvector was selected as the initial guess for a new calculation 
thereby speeding the computation. The output of the MATSUB program was 
also modified so that the eigenvectors were arranged into two nxn 
matrices corresponding to the real and imaginary parts of the transfor- 


mation matrix G, 
18 








From the two n x n matrices the 2n x én matrix H was formed and the 
inverse taken. The canonical state variables were then determined and 
the proper control calculated and applied. 

The path of the system trajectory in the physical space (x) was 
then calculated for a selected time interval between sampling instants 
by the Runge-Kutta-Gill method. Instantaneous values of the elements 
of the coefficient matrix F 2 F(x) were calculated for the new sam-= 
pling point and the procedure repeated until the trajectory in the 
canonical space was within a prexset epsilon region of the origin or 
a given time had elapsed. 

Initially the ideal case was simulated by holding the solution of 
the system trajectory at the sampling point until the new calculation 
of the control function was made and the control applied. Simulation of 
the non-ideal case, a finite calculation time, was affected by applying 
the control calculated at the (n-1)®" sampling, at the time of the nth 
sampling. This is equivalent to a controller continually applying a 
new control as soon as it can perform the sampling and control cal- 
culations. 

Simplified flow charts of both the ideal and non-ideal case simula- 
tions are presented in Appendix I. The computer program in FORTRAN 


language (for the ideal case) is also included in Appendix I. 
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4, Simulation of Control of a Non-Linear Second Order System. 
A. Background. 
The non-linear VanderPol equation 

% 2) 2x2) oar wo (4-1) 
was selected as the system on which to evaluate the controller. Al- 
though it is only a second order equation it was felt that since it 
is quite non-linear it would be a fair test of the controller and the 
two dimension phase space has the advantage of making graphical analysis 
more simple. 


The equation exhibits a stable limit cycle of the general shape 


shown in Fig. (4-1). 


Fig. (4-1) Phase Portrait For an Uncontrolled VanderPol Equation. 

















Referring to (4-1), when x = O, the system is most unstable. For 
this point the roots are +.5 + j866. It was determined from trajec~ 
tories of the uncontrolled system obtained by digital computer analysis 
that when the system was in the limit cycle the maximum excursion of 
displacement was ~~ 2.1925. This value was substantiated by an analog 
computer simulation. From (4-1) the eigenvalues, or roots, of the system 
when the displacement reaches its maximum value in the limit cycle were 
calculated to be -3.5235 and -.2835 so that the locus of roots of equa- 
tion (4-1) during a complete circuit of its stable limit cycle is shown 


in Fig. (4-3). ice 


5 +},.866 


~3.5235 ~e 2335 +o 
-3 -2 “4 | 


Sn-4 866 


Fig. (4-3) Root Locus of the Uncontrolled VanderPol Equation 
in Stable Limit Cycle. 

Examination of Fig. (4-1) or equation (4-1) reveals that the system be- 
comes unstable whenever |x| < 1.0 and becomes stable again whenever 
|x|> 1.0. Thus the equation presents the controller with the problem 
of controlling a system which contains at various times, stable real 
roots, stable complex roots, and unstable complex rootse 

It is informative to observe that a stability analysis of the 
VanderPol equation by Lyanunov's Second Method yields precisely the 
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same results as the more familiar methods of analysis such as the root 
locus. 


Rewrite (4-1) in vector form. Then 


wo lies 

=F fa (1-x)) i (4-2) 
or 

xX, = Xp (4-3) 
and 

Xo = -Xy + (1ex?)x, e (4-4) 

Choose 

vx) = xf = x + % 
then 

V(x) = 2x)x) + 2XoXo (4-5) 


and from (4-3) and (44), 
Vix) = 2X1Xo + 2Xo(-x, + (Ex 

or 

V(x) = ox8 (1-x‘’) (4-6) 
so that V(x) is positive definite but -V(x) is not for |x,|<1. There- 
fore the region near the origin (-1<x<1) is unstable. This is pre- 
cisely the result obtained by observing that the root locus passes into 
the right half plane whenever |x| <1. 

B. Procedure. 

It was decided to initially investigate the action of the control 
with no time delay between calculation of the control and its application 
in order to verify that the method of control and the simulation pro- 
cedure (described in Section 3.) were feasible. First a family of tra- 
jectories of the system with no control effort was obtained from 
various initial conditions chosen inside and outside the limit cycle 
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(Graphs A-1l through A-8). Trajectories were then obtained using the same 
initial conditions with the control system activated (Graphs B-1 through 
B-8). The sampling rate was set at e3 seconds and u was set at 1.0. An 
arbitrary epsilon region surrounding the origin was established so that 
the trajectory was considered to have reached the origin when the sum 
of the canonical state variables squared was less than .001. Whenever 
the trajectory entered the epsilon region of the origin the solution 
was terminated. 

Effects on the trajectories due to varying sampling rates (Graphs 
C-1 through C-5) were then investigated, and finally two trajectories were 
obtained for the system simulating the effect of finite calculation time 
(Graphs D-1 and D-2). These trajectories all were started from an 
initial condition close to the uncontrolled limit cycle trajectory. 

The simulation program included a sub-rottine for the graphical pre- 
sentation of the system trajectories. The graphs noted above are pre- 


sented in Appendix II. 


24 











ud 


ul 


a 


ay 


th 


as 


°norczed pucass guatq ‘septose 
GPWIT SL4ASTIEQZCCAIBYO peaz)equry 





S40 Vite 


nde 





Rn ee a ee 


tb 


ud 


bh 


a 


Cc 6 
¢ 99 
tro “tr 


& 
ca CS 


Sia. 5 gramainarenatao  w 





¢ 
an °S 


SEATVANGD Ta 
ea vis 


LSV2T 





Lot. 
GL8° 17 


¥S0"~ 
¥6°ot~ 


7 a 
TO8°s~ 


o2b6°™ 
O39" 2" 


inge* = 
S49°o~ 


SGn0-% 


59° s> 


Siege 
509° S" 


SlT¥A 
> “Neola 
< trys 

ELSON 


Sse oe 


CKOOaS £°O = 














O%o™ 
OF S= 
> 
= 
m3 
Qo: 
— O°? 
ry 
iy 
te 
{4 
a2 0°¢ 
o3 
fe) 
O72 
o°t 
;* 
a! 
NITOYEO ne 
HOFER 
OL SHOT LIGii100 
SCNOORS IVI ZINI 
GLVEY INITIALS : 








°SShTUSSY UCTABTNETS LO 


GV 


L7¥ 


g-¥ 


Gv¥ 


vy 





‘AS¥O GeTIOdLNOONN “TI 








eo 


UOTIG [RGU] 


79) 


25 














tS 










a ui eS ae ~ come Od Gag= 
Us ig 
3 hl. C2a = mais O°2¢ Oss 
6S0°> 
7 %6°9T° 99°2eT 0°? 0°% 
“ uu ee 0°? O°Ss 
uo $3°OT OFZ O°¢ 
rf $9 °OT G°T G°T 
GL96°0 F 
$3 F9LS° 82° Lt G° Ge 
gog°f = g* §988° ‘. 
xeage peddogs uoTantos 9ag°f + ge ¥9% as a as 
SUNTVANGO TM” ot] « SER FRA NIDINO i x YaGinn 
SME WIGN TICVLS “Nad HOVHY t— Ndvuy 
ESVE CURR AS OL SHOILIGHGO 
ISON SGNOODES TVI LINT 





aes a PETE agp Teal 




















—*O°T = un ‘"Squogss O° = GLVE ONITGYYS ‘'SV¥O GSTIOLLNCO °VeIT 


26 


















Cue 














ti g19%e~ 92°OT 
Led a @ = 
7 GLeL° 69°32 
u 6y° LT 
2 
- 269°9> S7°Q 
G0S¢’> 
v CO Rt 
a i ee Od 
tope°t - g6egse 
" twyerf + segee° 9°9 
GTULT GUTG ROTANTOS UTUQTM UOTeed g,ie°f - 93eP° 
Us{Tsds pexzequo’ Soysogcer easy TLV GL4LB°F + Bsr? eat 
con! on SEQUTVANSOTE |. | P -SGNTVANEDIa NIDINO 
SyUviay PIAVLS TIVLS HOVaAY 
ISVXI LS OF OL SNOILIGNOO 
canoogs IVI LINI 

















*GAL = on ‘siveguses7e 


at 
fx 
te 
az 
ayer 


L ONIINVS %SSVO GaTiOWINOD *a-II 


ot 








‘a*y WOT_oeg Bog *THOoUDIIE ST UOTANTOS STUL 











teeye~ 
! — 
alt ; go°g 








*spucoces W904 
“JtJ z0o1yzv poddoys uotgntog 


ybLa*- 
S3p?t- 
9rag°- 
ees a 

















santyaxcota | SGN TVAN@OTa HI DINO 

SRT TIEVLS TIEVIS HOVEY 
ISVE1 ESO OL 

SaNOOgS 











> eer we 














“O°T = mR "2°*t = ¥ 3°T = SNOTLIGHOO IVILINI ‘*2S¥O daTiOW¥LNod 


= 


Ga, 


YaGION 
Hava 





a aS 


28 











"oToAO 7 

GTTITT POTIJTpOW B sr1EeqUT 4 
perp os oe 
eog7: @ 3s 





wHaOn 























NIOIYO (SCNOOGS ) 
Suavring TL1EVLS TIGVLS HOVGY qiVe H2 yu 
LSVaTI ESOH ONT TERYS 








“OPT oe NM fG*T SX 6fG*T FX 6 SHOTLIGNOO IVICINT ‘asvO AYIGG"SL PAT 

















ra, 








D. Discussion of Digital Computer Simulation. 

Comparison of the trajectories of the controlled and uncontrolled 
systems starting from the same initial conditions (Graphs A-1 through 
A-8 and B-1] through B-8) clearly indicates that the control system with 
a control effort = 1.0 and a sampling rate of .350 secondsit is capable 
of driving the system into the unstable region surrounding the origin and 
maintaining it within some region of lesser extent than that enclosed by 
the system's uncontrolled limit cycle. Referring to Fig. (1-1), if the 
radius of S(R) were established for this system as R = .5, then the con- 
trolled system might be defined as stable while the uncontrolled system 
would be defined as unstable. Graph B-l, initial conditions x = .2, 

X = .2, illustrates the system trajectory in the region close to the 
origin. The chatter-motion mode of operation is readily apparent. 

The theoretical analysis of Section 2., predicting that by employ- 
ment of this type of controller, an unstable system could be made stable, 
but not asymptotically stable, and that a chatter-motion mode of 
operation would result around the origin, is substantiated. 

Examination of Graphs C-l through C-5 reveals that an increase in the 
control system sampling rate causes the trajectory, having once arrived in 
the neighborhood of the. origin, to remain thereafter within a smaller re- 


gion of the origin than the same system operating with a low sampling rate, 


Curves B-3 and C-3, both trajectories for the identical initial condi- 


tions and system parameters, demonstrate remarkably different trajectories, 


each of which eventually arrives in a neighborhood close to the origin. 


#The time interval between successive calculations of the control 
function. The Runge-Kutta integration of the trajectory employed an 
interval of 0.03 seconds. 
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Examination of the calculated data revealed that in the case of 
Graph B-3 the initial values of the canonical state variables were 
Y, = 075 = jol561 
Yo = 75 4 je1561 
while in the case of Graph C-3 the values were 
Y, = 0f5 =5-1561 
Yo = ¢75 <-je1561 . 

Since the control function u = -1.0 sgn ( DRy; + cd Y3)>s 
it is evident that u assumed a different sign in each case. That this 
was so can be determined by examination of the diverging paths of the 
trajectories in the two cases. 

It seemed evident that the more direct trajectory should be the 
correct one, however the possibility was raised that there might not 
be a unique solution for the control function. Six additional runs 
each exactly coincided with the more direct trajectory B-3, indicating 
that verhaps an error in calculation led to the selection of the op- 
posite value of control function rather than there being two possible 
solutions for the control function. 

The computation routine had been programmed so that the inverse 
transformation matrix was printed at each sampling point. It was dis- 
covered that the transformation matrices in the two cases were 


not identical. 


INITIAL INVERSE TRANSFORMATION MATRIX FOR GRAPH B-3 


o30820B-10 ©54051E-00 e S0000E=00 ¢ 400528 -00 
- 64051E-00 055820E-10 -.40032E-00 .50000E~00 
= 250820E-10 -=.64051E-00 .50000E-00 -.40032E-00 

0640516-00 -.35820E-10 °40032B-00 e SOO000E=00 
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INITIAL INVERSE TRANSFORMATION MATRIX FOR GRAPH Cx3 
e57180E-08 .64051E-00 .50000E-00 .40032E=-00 
=e64051E-00 .37180E=-08 -.40032E-00 .50000E5-00 
e50000E-00 .40032E-00 .65670E-09 .64051%=00 
~240032E-00 .50000E-00 -.64051E-00 .65484E.~09 
In order to determine if two solutions were possible, the two in- 
verse matrices were checked against the original matrix to see if both 
were valid. 
The original matrix was calculated by hand in the same manner that 


the computer was programmed to doe 


The procedure is illustrated below. 


The initial F matrix was lo, ; a and the initial eigenvalues were 


~.625 + j.7806. Then for the first eigenvalue 


0 u x - ; x 
EK 11 | BH a (7625 + 57806) Hi 


or 
X, ae (605 jo 7806 )x, 
x) = (+.625 - j-7806)x, 
so that 
Xo = La O 
Xx} = =socn =s jo 7806 © 


In a similar manner, for the second sigenvalue 
Xo od slay, 
X, = ~0625 4 je/806 . 


Then the 2n x 2n matrix becomes 


=0625  -.7806 ~.625 - 7806 
©7806 ~=.625  =.7806 ~.625 
1 0 1 0 
0 1 0 1 fo 


When this matrix is multiplied by the inverses, B-3 and C~3, only the 
inverse B-3 results in the identity matrix, demonstrating that the 
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calculation resulting in C~3 is not an independent correct solution. 
Comparison of the false trajectory with the correct one indicated 
the possibility of there being at least two other instances of erroneous 
calculations besides the one et the origin. These are indicated by the 
abrupt changes in direction (outward) of the trajectory at the points 
x = 65, y = 0.6, andx s <.5, y = .45 . Calculations similar 
to the one carried out above verified the fact that erroneous determina~ 
tions of the correct control function had been made. 
The similarity of the trajectories of curves C-l1, Ce2 and D-1l in 
the fourth quadrant raises the interesting comparison between the control 
obtained by Lyapunov's method and the method of optimum control employing 
the calculation of the optimum switching lines in negative time as dis- 
cussed by I. Flugge-Lotz and HA. Titus’. 
We may speculate that there is some optimum switching line for 


this system, emanating from the origin in a manner similar to that 


in Fig. (4-4). x 


Fig. (4-4) System Optimum Switching Line. 








The optimun trajectory for this system would then proceed from the 
initial condition to its first intersecticn with the optimum switching 
line, change the sign of its controller, and proceed along the optimum 
switching line to the origin. 

It appears evident that the method of switching logic discussed in 
this paper selects a switching time too soon in the cases of curves C1 
and C-2. In the case of curve De2 the system appeared to be in a state 
of indecision. The switchings at points A and C undoubtedly would have 
brought the trajectory initially closer to the origin than the one 
finally selected at BE. The switchings at B and D served only to hinder 
the system response. Investigation of the computer data again indicated 
erroneous calculations made for switchings at points B, C and D, although 
the resulting decision at C proved to be a correct one. 

Graphs D-1 and D-2, showing the simulation of the system operating 
under non-ideal conditions, again substantiate the fact that higher 
sampling rates generally produce more satisfactory operation of the 
system. It is interesting to note that the controlled system operating 
with a time delay of .3 seconds eventually settles into a stable limit 
cycle of approximately one half the size of the uncontrolled system. 

The effect of the controller on the excursions of the system roots 
as the state point travels along a trajectory for a particular initial 
condition is demonstrated in Figs. (4-5) and (4-6). The numbers indicate 
time in seconds to reach the position indicated. 

The root locus for the uncontrolled system starting from initial 
conditions x = .5, x = .5 is illustrated in Fig. (4-5). The roots 
attain their most stable positions after the trajectory has settled into 


its stable limit cycle. 
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Fige (4-5) Root Locus for Uncontrolled System, Initial 

Conditions x = oD, x = ole 
As the controlled trajectory approaches the region of chatter- 
motion surrounding the origin, the root locations depart very little 
from their most unstable position (4.5 + j.866) because maximum 

instability occurs when x = OO, The most stable root positions on 


the controlled system root locus occur early in the solution while 


the trajectory is relatively far from the origin. 


Fig. (4-6) Root Locus for Controlled System, Initial 
Conditions = 2925, x & =o. 
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Se Conclusions. 

The operation of the controller in regulating a simulated non- 
linear system described by the VanderPol equation supports the theo- 
retical predictions developed in the design. The solution obtained by 
this method is seen to be less than optimum because tho logic does not 
cause the control to switch at the optimum switching curve. It is 
apparent from the figures that the sampling rate for determining the 
eigenvalues, and hence changes in control action, is critical for this 
to be an effective control design procedure. With the fine sampling 
the system disturbances were rapidly and effectively controlled. 

There exist today very few direct techniques for the control of 
non-linear systems. The procedures studied here provide such a tech- 
nique. The method has the advantage of being applicable to non- 
linear time varying systems of any order. It does however require a 


digital computer in the control link. 
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APPENDIX I 


FLOW CHARTS AND FORTRAN LANGUAGE PROGRAM 
FOR DIGITAL COMPUTER SIMULATION 


of 
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Fige I=-1 FLOW CHART FOR IDEAL CASE SIMULATION. 
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Fig. I-2 FLOW CHART FOR FINITE CALCULATION TIME SIMULATION 
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